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ABSTRACT 

Context. In the context of planet formation, anticyclonic vortices have recently received lots of attention for the role they can play in 
planetesimals formation. Radial migration of intermediate size solids toward the central star may prevent their growth to larger solid 
grains. On the other hand, vortices can trap the dust and accelerate this growth, counteracting fast radial transport. Multiple effects 
have been shown to affect this scenario, such as vortex migration or decay. 

Aims. The aim of this paper is to study the formation of vortices by the Rossby wave instability and their long term evolution in a full 
three dimensional protoplanetary disc. 

Methods. We use a robust numerical scheme combined with adaptive mesh refinement in cylindrical coordinates, allowing to afford- 
ably compute long term 3D evolutions. We consider a full disc stratified both radially and vertically that is prone to formation of 
vortices by the Rossby wave instability. 

Results. We show that the 3D Rossby vortices grow and survive over hundreds of years without migration. The localized overdensity 
which initiated the instability and vortex formation survives the growth of the Rossby wave instability for very long times. When the 
vortices are no longer sustained by the Rossby wave instability, their shape changes toward more elliptical vortices. This allows them 
to survive shear-driven destruction, but they may be prone to elliptical instability and slow decay. 

Conclusions. When the conditions for growing Rossby wave-related instabilities are maintained in the disc, large-scale vortices can 
survive over very long timescales and may be able to concentrate solids. 

Key words. Planets and satellites: formation - Protoplanetary discs - Hydrodynamics - Instabilities - Accretion discs 



1. Introduction 

The origin of the kilometer size planetesimals is still an issue 
of planet formation theory. The growth of micron size solids 
toward centimetre or meter blocks is predicted by coagulation 
models (Dominik 2009 ), but collisions usually lead to destruc- 
tion of such solids (Benz 2000). These collisions arise because 
the pressure supported gas moves at sub-Keplerian speed and 
the solid bodies experience a head wind from the gas due to the 
drag forces. The consequence is a loss of angular momentum 
and a radial drift of solids spiraling toward the central star. This 
can occur on a timescale as short as a hundred years for meter 
size blocks at one astronomical unit ( fWeidensc hilling 1977]). 
Such timescales are far too short to explain subsequent growth 
into larger particles that are unaffected by the head wind. 

To overcome these difficulties, anticyclonic vortices, where 
the velocity streamlines rotate in the opposite direction to the 
1990| ), may be a good environment 



Keplerian flow (jMarcus 



for planetesimal growth. They induce a drag force towards the 
centre of the structures, and are therefore proposed as possible 
nurseries for intermediate size blocks, concentrating the solids 
and accelerating the planetesimal formation processes (Barge & 



Sommeria||1995 


Tanga et al.||1996l |Bracco et al.||1999[|Godon 


& Livio 1999 


2000, Johansen et al. 2004 Heng & Kenyon 


2010). As mentioned by |Armitage| (|201 1\, vortices can also 



be important for mechanisms that require an increase in the 
dust-to-gas ratio, such as the streaming instabil ity (|Johansen| 
|et al.|2009] ). Furthermore, |Klahr & Bod enheimer ( 2006]) argued 
that anticyclonic vortices are regions with lower turbulence. 
Fragmentation is less likely to occur if lower velocity fluctua- 
tions prevail. 

One of the main difficulties with this scenario is the question 
of the formation of vortices. Different instabilities have been 
proposed for the vortex generation, such as the baroclinic 
instability ( |Klahr & Bodenheimer|[2003] |Hahr||2004|) in which 



interest has been recently revived ( |Lesur & Papaloizou |2010 



Lyra & Klahr 201 1\, or poten tially the magneto-rotational 



instability (Fromang & Nelson] [2005). In this paper we will 
study the Rossby wave instability (RWI), that may form vortices 
with unusual elongated shape in regions of particular interest for 
the long term survival of those structures. It has been recently 
pointed out that the survival of the vortices on a sufficient 
timescale may be an issue. The elliptical instability may destroy 
three dimensional (3D) elliptical vortic es (|Lesur & T apaloizou 
2009). Moreover, Paardekooper et al. (2010) have shown that 



vortices can be subject to radial migration toward the disc centre 
in a short timescale. 

In this paper, we investigate the generation, the 3D struc- 
ture and the long term evolution of vortices formed by the RWI 
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Fig. 1. Mid-plane density of the gas in g cm 3 . The density bump 
is placed at 3 AU. 



within a protoplanetary disc. This follows the work presented in 
Meh eut et al.| ( |2010| l, revisited with a new code that allows to 
handle long term simulations. The numerical and physical setup 
is described in the next section. The results are presented in sec- 
tion^ followed by a detailed discussion including comparisons 
with previous works. 



2. Methods and setup 

2. 1 . Rossby wave instability 



The RWI ( jLovelace et al.|[l999l |Li et al.|2000] [2001 ) has been 
studied and discussed in various situations of differentially ro- 
tating discs, from gala ctic discs ( |Lovelace & H ohlfeld 1978; 
Sellwood & Kahn 1991]> to m icroquasar and protoplanetary discs 
(Papaloizou & Pringle 1985 ; Lovelace et al.|1999] >. It can be seen 
as the form that the Kelvin-Helmholtz instability takes in differ- 
entially rotating disks, and has a similar instability criterion: It 
requires an extremum in a vorticity related quantity, defined in a 
non-magnetised thin disc as ( [Li et al.|2000| >: 



(pE-r)2/r 



2(V x v). 



(1) 



where S is the surface density, p the pressure, v is the velocity 
of the fluid, y the adiabatic index, £2 the rotation frequency and 
k 2 = 4£2 2 +2r£2£2' the squared epicyclic frequency (so that k 2 /2Q 
is the vorticity). Here the prime denotes a radial derivative. For 
the isentropic discs we consider here, this criterium is reduced to 
an extremum of . This quantity is directly related to vortensity 
defined as the ratio of vorticity to surface density. 



2.2. Equations 

We work in cylindrical coordinates (r, z, f) with the 3D Euler 
equations 



d,p + V • (vp) = , 
d,(pv) + V ■ (vpv) + Vp = -pVO G , 



(2) 
(3) 



where p is the mass density of the fluid, v its velocity, and p its 



pressure. 



ton = - 



CM., 



is the gravity potential of the central 
object with G the gravitational constant and M, the mass of the 
star. We consider a barotropic flow, i.e. the entropy is constant in 



the entire system. The pressure is then p = Sp 7 , with the adia- 
batic index y — 5/3 and the constant S related to entropy. The 
sound speed is given by c 2 s = yp/p = Syp y ~ l and the tempera- 
ture by T ~ pip — Sp y ~ l . The temperature is normalised to the 
temperature at 1 AU. 



2.3. Numerics 



The numerical methods used for these simulations are inspired 
by those of M eheut et al.| ( |2010) . We use the Message Passing 
Interface -Adaptive Mesh Refinement Versatile Advection Code 
(MPI-AMR VAC) developed by Keppens and Meliani ( |Keppens 
et al.|2012) . The use of a code that allows the mesh to adapt dur- 
ing the simulation has multiple advantages. One aim is to reach 
higher resolution in the vortex regions, as well as to have a better 
resolution in the upper region of the disc, where the gas density 
decreases abruptly forming a 'corona' . The use of AMR offers 
flexibility to enforce a lower resolution in this physically irrel- 
evant, but computationally challenging, region. Since low den- 
sity material can easily be accelerated, in turn enforcing smaller 
timesteps to maintain numerical stability in explicit time step- 
ping schemes, handling the corona at low resolution is computa- 
tionally advantageous. The current AMR approach can therefore 
model the unstable disc on long timescales with similar comput- 
ing resources as the short timescale run done earlier. Also, by 
better numerically representing the disc-corona interface, com- 
bined with a sharper limiter function, we have a more robust 
numerical treatment of the governing equations. 



The numerical scheme is the same for all refinement levels, 
namely the Total Variation Diminishing Lax-Friedrich scheme 
(see Tdth & Odstrcil (1996 )) with a third order accurate Koren 
limiter ( |Koren|199'3] ) on the primitive variables. We use a cylin- 
drical grid with r e [1,6] AU, z e [0,0.5] AU and the full az- 
imuthal direction <p e [0, 2n]. The simulations consider only the 
upper half of the disc, as the disc mid-plane appeared to be a 
symmetry plane for the instability in our earlier full vertical sim- 
ulations (Meheut et al. 2010[). If not specified, the length is given 



in AU and the code unit time corresponds to l/(2n) yr. The res- 
olution of the base level is (64, 32, 32). Up to three levels of re- 
finement are allowed, corresponding to an effective resolution of 
(256, 128, 128). In the region of the density bump where the in- 
stability is expected to grow, the resolution is fixed to the higher 
level during the whole simulation. The initial grid is presented 
in Fig. [2] and then evolves to follow the growth of the instabil- 
ity. The refinement criterion being the density variation, the grid 
follows the growth of the spiral density waves propagating radi- 
ally. At the end of the simulation ~ 70% of the grid volume is 
at the highest resolution level. In the beginning, it is only about 
~ 26%, as visually clear from Fig. [2] Since refinement happens 
in all directions, each grid block handled at the base resolution 
level represents a gain factor of 64 in computational cost com- 
pared to a uniform high resolution run. This gain is even more 
dramatic when memory issues are considered as well. During the 
entire run, the corona region is maintained at this lowest resolu- 
tion, dramatically impacting stability, memory and wall-clock 
time. 

The boundary conditions are transparent at inner and outer 
radius, we consider a mid-plane symmetry and transparent 
boundary conditions are also applied at the upper boundary, 
whereas the azimuthal direction is periodic. 
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Fig. 2. Initial gas density in the disc in gcrrT 3 (up) and AMR grid (down). The three levels are plotted in white, blue and dark blue 
from lower to higher resolution. The highest resolution is reached in the overdensity region and at the interface between the disc 
and corona. 



2.4. Initial conditions 

The initial conditions for the simulations are chosen to represent 
a protoplanetary disc in radial and vertical equilibrium. We con- 
sider that an overdensity is formed at some radius (rbump) of the 
disc. This bump could have been formed due to multiple effects. 



1. 



The presence of a dead zone in the protoplanetary disc, cor- 
responding to a region of lower ionisation and resistivity, 
leads to a lower accretion rate in this region (Gammie 1996). 
This can induce the formation of an overdensity at the edges 
of this region (|Varniere & Ta gger 2006} |Lyra et aL |2008 



Kretke et al.|2009) >, which can in turn trigger the RWI. The 



initial density bump used in this paper, is a Gaussian fit to 
the one obtained in Varniere & Tagger (2006), when the in- 



stability started to grow. 
The ice sublimation front ('snow line') can also be responsi- 
ble for the formation of an overdensity (Kret ke & Lin|2007] l. 
One can note that an extremum of entropy should be located 
in this region, which may also trigger the RWI (Lovelace 
|et al.|1999> . 

It has also been shown that the RWI grows at the edge of 
planet gaps ( Koller et al.|2003 Lin & Papaloizou|2011| l. 



The radial positions of these regions depend on the physical 
characteristics of the disc (such as accretion rate and tempera- 
ture) or the position of the planets. We have chosen to place the 
bump at a distance of 3 AU which is a plausible region for these 
different effects. 

The initial mid-plane density is given by 



(4) 



gas is then sub-Keplerian. There is no vertical velocity initially, 
whereas a very low radial velocity perturbation is added to the 
equilibrium as a seed for the instability. These are random per- 
turbations, meaning that all the spatial frequencies in the three 
directions are present in the seed. 

3. Results 

The evolution of the pressure bump in a protoplanetary disc is 
followed over more than 600 years (t ma x = 4000). This timescale 
corresponds to ~ 123 rotations of the overdensity and ~ 43 ro- 
tations at the outer edge of the grid. Such a configuration allows 
the development of the Rossby wave instability and we present 
here the characteristics of the 3D flow under this instability and 
its non-linear evolution. 

3.1. Growth of the RWI 

The RWI is characterised by vortex waves in the region of the 
pressure bump and spiral density waves propagating on each 
side of the bump. These features can be seen on Fig. [3] where 
the vortex waves can be identified by the vortensity extrema or 
the velocity streamlines characteristic for such waves. On the 
upper plot, the inner and outer Rossby waves are visible. Spiral 
density waves are emitted from each vortex inward and outward 
(e.g. at t — 300). The growth of the instability is quantified in 
Fig. [4] and [5] The linear phase of the instability corresponds to 
exponential growth of the perturbations until t ~ 300 when the 
density perturbation reaches ~ 20% of the initial density. The 
straight line corresponds to a linear fit of the amplitude growth 
of the density perturbation in logarithmic scale. This fit gives a 
growth rate of y — 0.033Qo in code units which corresponds to 
0.17Q(r B ) or 0.21 yr _1 . We have also plotted the amplitude evo- 
lution of selected azimuthal modes, out of the Fourier transform 
of the density 



with po = 3.10~ y gcirT 3 the density at r — 1 AU, and a = -1.5 
the power law index of the underlying density. This value for the 
alpha parameter gives a surface density varying approximately 

as I ~ r l/2 in the absence of the bump. Parameters r B = 3 AU, P( r ' z ' <P> f ) = 2j Pm(r ' z ' f) ex P( _I ' m V)- 
X = 1, <t = 0.1 AU are the radius, amplitude and width of the 
Gaussian bump. 

The vertical and radial force equilibria give the initial verti- 
cal structure in density and azimuthal velocity, respectively. The 
pressure is calculated with constant entropy, with S - 10~ 3 . The 
mid-plane and vertical density profile chosen as initial condition 
are shown on Fig. [T] and Fig. [2] 

The azimuthal velocity is close to Keplerian where the de- 
viation from Keplerian rotation is due to the pressure gradient. 
In most of the disc, except in the inner part of the bump, the 



(5) 



The most unstable mode during the exponential growth has az- 
imuthal mode number m = 5, but lower azimuthal mode num- 
bers are then excited by the m = 5 mode leading to their growth, 
with m = 1 and m = 2 shown on the figure. The mode number 
corresponds, on Fig. [3] to the number of anticyclonic vortices, 
which are those counter-rotating the Keplerian rotation. They 
correspond to high pressure and negative vortensity regions. A 
high number of anticyclonic vortices come out from the initial 
random perturbations, rapidly decreasing to 5. 
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Fig. 3. Density (left), perturbation of density and velocity streamlines (center), perturbed vortensity defined as £ - ) tp with ^ 
(V x v)Jp (right) at t = 100, 300, 800 from top to bottom. 



3.2. Saturation of the RWI 

After the exponential growth, the instability reaches saturation 
due to non-linearities inducing mode mixing. Fig. [6] shows the 
amplitude of the modes at saturation. The fifth mode clearly 
dominates for quite some time, and modes with high azimuthal 
mode number are less prominent, but very low and intermediate 
azimuthal mode numbers are non-linearly excited by the prevail- 
ing mode. However one can see on Fig. |3]that later on (t = 800) 
m = 2 dominates, and then m = 1 until the end of the sim- 
ulation (see Fig. [6]). On the other hand, the higher mode num- 
bers (5 < m < 10) are not dominant after saturation. So the 
general sketch of the instability is first the linear growth of the 
fifth mode, then saturation is reached and the larger wavelength 
modes (lower mode number) dominate with a transfer towards 
larger scales. This evolution diminishes the number of vortices 
through merging as described in Godon & Livio|(|1999[ ). This is 
similar to the 2D simulations of Baty et al.| ( |2003| ) "where pla- 
nar shear layers were studied in high resolution Cartesian set- 
tings. They found that the pairing/merging of vortices in Kelvin- 
Helmholtz unstable shear layers, studied from hydro to magne- 
tized cases, is controlled by growth of subharmonic modes and 
phase differences between them. In the present 3D cylindrical 
configurations, their results represent local box evolutions, per- 



formed in the frame rotating at the local Keplerian speed. This 
pairing/merging trend to large scale structure is then predomi- 
nantly a 2D process known from planar hydrodynamical evo- 
lutions. As the RWI is closely related to the Kelvin-Helmholtz 
instability, the merging of vortices seen in our 3D, thin disc sim- 
ulations is consistent with the Baty et al. (2003 i findings. 



3.3. Decay of the vortices 

The simulation has been run over more than 600 years to inves- 
tigate the long term evolution of the vortices. After the growing 
phase, the vortices survive during hundreds of years but then 
start to decrease and have disappeared at the end of the simula- 
tion. After a few hundreds of inner rotations, the dominant az- 
imuthal mode does not change anymore, and has a global m = 1 
value. This corresponds to the largest scale possible. One can 
see on Fig. [7] that during the decaying phase (after about 100- 
150 years), vertical mode structure appears in the m = 5 mode. 
This structure is present inside the vortex (r - 3) but not in the 
outer region (r = 3.5). 

One can notice that the spiral density waves disappear af- 
ter the growing phase when the vortices survive. This implies 
that Rossby wave instability is not active anymore due to the 
weakening of the bump by the instability itself and no angular 
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Fig. 4. Time evolution of the amplitudes of the density pertur- 
bations in a logarithmic scale. The time is given in code units 
(lower axis) and in years (upper axis). We also show the am- 
plitude of some of the modes with low mode number. The upper 
figure is a zoom on the exponential growth where we have added 
a linear fit for the dominant mode whose growth rate is 0.033 in 
code units. 



Fig. 5. Time evolution of enstrophy in a logarithmic scale. The 
time is given in code units (lower axis) and in years (upper axis). 



4000 




Fig. 6. Comparison of the amplitude of the first 10 azimuthal 
modes between t — and t - 400 (left) and over the whole 
simulation (right). The instability is dominated by the mode m = 
5 but both lower and higher azimuthal modes are excited during 
the growth phase. The decaying phase shows a shift to lower 
azimuthal mode number. 



momentum is transferred radially. This is quantified in Fig. [8] 
where the radial accretion rate is plotted: the bump stops to de- 
crease due to radial accretion. When the Rossby vortices are not 
anymore sustained by the RWI, they should be very efficiently 
destroyed by differential rotation (Tagger 2001), but the shape 
of the streamlines changes to become more similar to closed el- 
liptical streamlines as presented in Fig. [9] The vortex streamlines 
are not directly linked to the shearing sub-Keplerian streamlines, 
and can survive over more than one rotation time. By the end 
of the simulation, the vortices have been completely destroyed. 
This will be discussed in the following section. 



4. Discussion 

4.1. Comparison with previous 3D simulations 

From a numerical point of view, the main difference with the pre- 
vious simulation of Me heut et al.| ( |2010| l is the use of an AMR 
grid. This allowed to increase the resolution at the interface be- 
tween the disk and the corona above it and to avoid numerical 
difficulties due to dynamics in this low density region. Our cur- 
rent simulation uses an improved limiter and benefits from grid 
adaptivity, allowing to robustly compute and follow the growth 
of the instability over longer times. The limit in time is now re- 
lated to computational resources. We tested simulations of the 
initial growth of the RWI with and without AMR, getting near- 
identical results. 
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Fig. 7. Vertical structure in density of the dominant fifth az- 
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Fig. 8. Accretion as a function of radius during the exponential 
growth {solid line) and at the end of the simulation {dashed line). 
The accretion rate is defined as - J J dipdzrpv r , and the ampli- 
tude is given in code units which correspond to ~ 10~ 4 M Q yr~ l . 





Fig. 9. Vortex streamline at t = 300 and t = 1000: the shape 
of the vortex evolves during the simulation. The radial extent 
of the vortices is constant and fixed by the radial extent of the 
density bump, whereas their azimuthal extent increases when the 
m mode number decreases. On the right figure, the streamlines 
become elliptic and an aspect ratio can be estimated (~ 6 here, 
but each vortex has a different aspect ratio as can be seen on Fig. 

0- 



4.2. Vortex migration 

We don't see any migration of the vortices in our simulation over 
more than 600 years. This result is expected, since the growth 
mechanism of these vortices is localised in the bump region. 
Vortex migration is oriented in the direction of positive density 
gradient, and as a consequence the presence of a density max- 
imum at the location of the vortices locks their radial position. 



One can see on Fig. 1 1 that the overdensity is diminished by the 
Rossby wave instability but it survives the growth of the RWI. 
The vortices are then indefinitely blocked in this region. 



The simulations presented in this paper show that the most 
unstable mode of the RWI in the current disc equilibrium has 
an azimuthal mode number m = 5. In Me heut et al.| ( |2010| ) we 
stressed the global m — 1 mode. There, only one mode was 
present in the seed perturbation, whereas we now choose to use 
random perturbations that include all modes. For confirmation, 
we have run the simulation of Mehe uFet al.| ((2010 > with ran- 
dom initial perturbations, and the instability was dominated by 
the m = 5 mode. The difference with our earlier work is then 
not due to disparity in initial equilibrium or numerical approach. 
This result is also coherent with linear calculations of Meheut 
|etaT1 < |20T2> . 

The unusual 3D structure of the Rossby vortices with an non- 
negligible vertical velocity is confirmed by these higher resolu- 
tion simulations, as can be seen on Fig. 10 We have selected 



some streamlines which pass next to the center of the vortices, 
they show the 'eye' of the vortices and the larger basis. The di- 
rection of the flow in the eye of the vortices differs between the 
cyclones (downward flow) and anticyclones (upward flow). In 
the outskirts of the vortices, the flows have opposite directions. 



4.3. Elliptical instability 

The elliptical instability is a purely 3D instability that is due to 
resonances between the vortex turnover frequency and an iner- 
tial wave frequency (Kerswell 2002). Lesur & Papaloizou (2009) 
have shown that 3D elliptical vortices are vulnerable to this in- 
stability, and can destroy them. In our simulation, the elliptical 
instability does not prevent the emergence of vortices by the 
RWI and their survival over hundreds of years. This was ex- 
pected as closed velocity streamlines are needed for the turnover 
frequency to be properly defined and the elliptical instability to 
grow. Some streamlines of the (non-steady) Rossby vortices can 
be seen on Fig.|9](see also Fig. 10 1. After saturation of the RWI, 
the streamlines tend to close, on the one hand this allows them 
to survive against shearing, but on the other hand, they become 
prone to the elliptical instability. This may be the reason for 
the decay of the vortices. The growth of this instability can ex- 
plain the appearance of smaller scale vertical perturbations (as 
the ones plotted in Pig. |7]» whereas the azimuthal scale is fixed. 
Indeed contrary to the RWI, the elliptical instability is a local in- 
stability without any interaction between the vortices, it can not 
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Fig. 10. Some streamlines which pass next to the centre of the 
anticyclone (upper figure) and cyclone (lower figure) at t — 300 
are plotted in 3D. The direction of the flow is given by the colors: 
the flow moves from the blue part to the red part. The cyclone is 
a downward flow, contrary to the upward anticyclonic flow. The 
vortices have a coherent structure over the whole height of the 
disc. The radial extent of this figure corresponds approximately 
to half of the one of Fig|9j only the inner part of the vortices are 
shown here. 



change their number but affects only their structure. The detailed 
study of the decay of the vortices by the elliptical instability, in- 
cluding elongated vortices or short wavelengths, needs a very 
high resolution as has been done in Lesur & Papaloizou (2009). 
They studied the effects of the aspect ratio, the ratio of semi- 
major to semi-minor axis of elliptical streamlines, and showed 
that elliptical vortices are always unstable. However, the ellip- 
tical instability becomes weak for large aspect ratios (i.e. elon- 
gated vortices, as is the case in Fig. [9}. Therefore, the low az- 
imuthal mode number vortices are expected to survive longer, 
as observed in our simulations 
|Livio 



As pointed out by Godon & 
( |1999| l, we must caution that the inherently higher numer 



Fig. 11. Mean azimuthal density in the mid-plane of the bump 
region, initially and at t — 1000. The bump has survived the 
growth of the RWI. 



buoyant, e.g. subject to a subcritical baroclinic instability, and 
Lesur & Papaloizou] d2010} have shown with higher resolution 
simulations that this secondary instability does not fully destroy 
the vortices (see also |Lyra & Klahr|201 1| >. 



ical dissipation of our finite volume method compared to a spec- 
tral method, also induces vortex decay. The elliptical instability 
is also present in the vortices formed when the disk is radially 



4.4. Outlooks 

We have shown that Rossby vortices can emerge in a 3D pro- 
toplanetary disc with an overdensity, and that they can survive 
for hundreds of years without migrating. When they are not sus- 
tained anymore by the RWI, the structure of the vortices will 
evolve to previously studied vortices in isolation ( Godon & Livio 
1999 Lesur & Papaloiz ou|2009| l. They are then expected to be 
destroyed and may be prone to the elliptical instability and vis- 
cous decay. However, our simulations show that the Rossby vor- 
tices can survive for long time scales, if they keep being sus- 
tained by the RWI and their unusual structure is preserved. That 
would be the case if the mechanism that forms the overdensity 
and launch the RWI, is permanent, such as present in a disc with 
both active and dead zone regions. 

A point to explore is the study of this instability when the ac- 
tive zone is also considered. On the one hand we would expect 
viscosity to decrease the growth rate of the instability and even- 
tually to destroy the vortices. On the other hand, if a dead zone 
is included, a density bump arises that will in turn sustain the 
instability. A full MHD simulation could allow to study jointly 
the bump formation process and its decay due to the RWI, in 
combination with magneto-rotational driven accretion. 

The results presented here are especially interesting in the 
context of planetesimal formation. Future work should also 
study the joint evolution of the gas and dust particles to ques- 
tion the influence of the dust particles on the growth of the 
RWI and their concentration in Rossby vortices. Previous studies 
have proposed a bi-dimensional appro ach (Chang & Oishi 2010; 
Johansen et al.|2004| |Wolf & Klahr|[2005| l but a full 3D study, 
as the one presented here including solid particles, will be nec- 
essary to handle both vertical stratification and vortex concen- 
tration. With this intent, a module of MPI-AMRVAC has been 
developed and was recently used for circumstellar wind model- 
ing ( |van Marie et al.|201 1) . 
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On an analytical point of view, the only work we are aware 
of, that deals with the RWI in 3D is the one by |UmurhaT?1 ( |20 1 0} . 



The vertical structure of the modes should be studied in more 
detail with dedicated tools as the one presented in|Zhang & Lai| 
([2006). This will be addressed in |Meheut etaL| ( |2012) . 

Acknowledgements. This work was partially supported by the Tournesol pro- 
gram of the PHC (Partenariat Hubert Curien) and by the Swiss National Science 
Foundation. 
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